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A  new  approach  for  the  formulation  of  boundary  conditions  for  the  counterflow  configu¬ 
ration  is  presented.  Upon  discretization  of  the  steady-state  Navier-Stokes  equations  at  the 
inflow  boundaries,  numerically  algebraic  equations  are  imposed  as  boundary  conditions, 
while  upon  discretization  of  the  unsteady  Navier-Stokes  equations  at  the  outflow,  differen¬ 
tial  boundaries  result.  It  is  demonstrated  that  the  resulting  numerical  differential-algebraic 
boundary  conditions  are  suitable  to  account  for  the  multi-directional  character  of  the  flow 
at  the  boundaries  of  the  counterflow  configuration. 


I.  Introduction 

The  counterfiow  configuration  provides  a  comprehensive  framework  for  studying  the  characteristics  of 
non-premixed  laminar  and  turbulent  flame  problems1- 1> .  However,  apart  from  the  simplified  one-dimensional 
spatial  models,  the  fidelity  of  direct  numerical  simulations  (DNSs)  for  the  counterflow  configuration  in  terms 
of  robustness/accuracy  exhibit  significant  sensitivity  to  the  boundary  condition  (BC)  treatment '  ~9 .  This 
behavior  is  mainly  due  to  the  multi-directional  character  of  the  flow  at  the  boundaries  which  must  be 
properly  accounted  by  the  BCs.  To  mitigate  this  problem,  Yoo  et  al.'  developed  improved  BCs  based  on 
the  Navier-Stokes  Characteristic  Boundary  Conditions10  (NSCBCs)  for  laminar  and  turbulent  counterflow 
flames.  The  major  improvement  consisted  in  introducing  the  (no-longer  negligible)  transverse  terms  into 
the  Locally  One  Dimensional  Inviscid  (LODI)  relations  in  order  to  capture  the  multi-dimensional  effects  at 
the  inflow/outflow  boundaries.  However,  two  major  shortcomings  can  also  be  identified  with  the  improved 
NSCBCs.  First,  they  preserve  the  use  of  relaxation  coefficients  entering  the  improved  LODI  relations  and 
these  coefficients  are  problem  specific.  These  coefficients  provide  an  optimal  balance  between  achieving  the 
prescribed  upstream  values  for  the  inflow  variables  and  reducing  spurious  wave  reflections,  and  usually  must 
be  determined  through  a  trial  and  error  process;  this  is  an  expensive,  time-consuming  procedure.  Second, 
once  derived  for  real-gas,  the  implementation  of  these  revised  NSCBCs  adds  a  considerable  computational 
cost. 

In  this  work,  we  adopt  a  totally  different  viewpoint  for  constructing  BCs  for  the  counterflow  configuration 
that  results  in  a  very  concise  framework.  In  particular,  we  combine  the  steady-state  Navier-Stokes  equations 
with  an  initial  flow  of  potential  type  to  construct  numerically  consistent  algebraic  BCs  at  the  inflow  bound¬ 
aries.  The  structure  of  the  paper  is  as  follows.  First,  the  governing  equations  of  the  problem  are  presented. 
A  detailed  discussion  on  the  construction  of  the  initial  profile  of  the  flow  follows.  Then,  the  new  BCs  are  pro¬ 
vided  for  the  inflow  and  outflow  boundaries.  The  specifics  regarding  the  numerical  implementation  that  was 
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employed  follows.  Next,  the  validation  of  the  proposed  differential-algebraic  BCs  with  numerical  examples  is 
performed.  Finally,  a  summary  is  given  and  several  aspects  regarding  the  proposed  approach  are  discussed. 

II.  Governing  equations 

A.  The  Navier-Stokes  equations 

The  Navier-Stokes  (NS)  equations  for  a  compressible  reacting  multicomponent  mixture  of  N  species  expressed 
in  terms  of  the  conservative  variables  U  =  [Um\  (m  =  1, . . . ,  N  +  4)  where 

U=  [p,pu1,pu2,pU3,pY1,...,pYN-i,pet]T  (1) 


and  in  Cartesian  coordinates  are 


dp  d(puj)  ^ 
dt  dxj 

d(puj)  d(pUjUj)  dp_  _  drjj 

dt  dxj  dxi  dxj 

d{pYn)  ,  d{pYnUj)  dJnj  ,  . 

dt  dxj  dxj  U’n 

d{pet)  d[(pet  +  p)uj }  =  _dqj_ 

dt  dxj  dxj 


(hi  =  1,2,3) 

(n  =  1, ...  ,1V  —  1) 

djTjjUi) 

dxj 


(2) 

(3) 

(4) 

(5) 


where  the  indices  i  and  j  follow  the  summation  convention,  t  is  time,  Xj  is  the  j-th  spatial  coordinate,  p 
is  the  mass  density  of  the  fluid,  Uj  is  the  j-th  component  of  the  velocity  vector,  Yn  is  the  mass  fraction  of 
species  n,  p  is  the  pressure,  et  =  e  +  \ui%ii  is  the  total  specific  energy  (e  stands  for  the  specific  internal 
energy),  is  the  j-th  component  of  the  mass  flux  vector  J„  of  species  n,  qj  is  the  j-th  component  of  the 
heat  flux  vector  q  and  uin  is  the  mass-production  rate  of  species  n.  Finally,  r;,  is  the  Newtonian  viscous 
stress  tensor 


/  dui  duj  2  duk  x  \ 
dxi  3  dxk  1J ) 


(6) 


where  p  is  the  dynamic  viscosity  of  the  mixture  and  6tj  is  the  Kronecker  delta.  In  the  present  equations, 
body  forces  have  been  neglected. 

Letting  F*  =  [F'm\  (m  =  l,...,A  +  4)  denote  the  flux  vector  of  the  conservative  variables  along  the  i-th 
direction,  i.e. 


F1  =  [pui,puiUi  +pSu,pu2Ui  +p62i,pu3Ui  +p53i,puiY1, . . . ,  puiYN-i,  (pet  +p)ui]T, 


(7) 


the  NS  equations  can  be  cast  in  compact  form  as 


<9U  dFl 
dt  dxi 


where  the  vector  C  =  [Cm\  (m  =  1, . . . ,  N  +  4)  is  defined  as 


C=  [0, 


dr-ij  dr2i  dr. 


2  j 


r3  j 


dJ 


1 3 


dJ, 


dx  ,■ 


dxj 


dxj 


dxj 


■  Wl,  ■ 


N-lj 


dxj 


+  UJn-I,  — 


dqj  d^ijUi) 


dxj 


dxn 


(8) 

(9) 


When  C  =  0,  one  obtains  the  compressible  Euler  equations  augmented  by  the  species  equations. 


B.  The  expressions  for  the  heat,  mass  fluxes  and  viscosity  of  the  mixture 

Here  we  employ  the  generalized  heat  and  mass  transport  equations  based  on  the  fluctuation-dissipation 
theory1 1-13 .  The  mixing  rules  employed  for  the  computation  of  the  mass  diffusion  coefficients  and  thermal 
diffusion  factors  have  been  derived  by  Harstad  &  Bellan14  and  are  based  on  the  coupling  of  nonequilibrium 
thermodynamics11  and  Grad’s  13-moment  theory15. 
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In  the  above  setting,  having  neglected  the  term  proportional  to  X7p/p  since  it  is  anticipated  that  its 
contribution  is  minimal,  the  mass  flux  vector  of  species  n  reads 


J  n  —  P 


Ad  X7T 

^(Auc  ~  +  YnDT,n^T 

.  k  ^ 


where  Mn  is  species  n  molar  mass,  Dnk  are  the  pairwise  mass  diffusion  coefficients  and 

Dnk  =  E  amkD%  DT,n  =  Y1  “T,kDkn, 


abT,r. 


-^fcaT,fen> 

k^n 


(10) 

(11) 

(12) 


where  Xk  is  the  mole  fraction  of  species  k,  a p  kn  are  the  binary  thermal  diffusion  factors  and  cionk  are  the 
mass  diffusion  factors  which  are  calculated  in  conjunction  with  the  EOS  as 


O'Dnk  — 


dXn 

dXk 


■X, 


dlQ7n 
;  dXk  ’ 


(13) 


with  7„  =  (t>n/(t>n,  where  <j>n  is  the  fugacity  coefficient  of  species  n  and  the  superscript  o  denotes  the  single 
species  (Xn  =  1)  limit.  The  elements  D^fn  are  obtained  as  the  solution  of  the  mixing  rules  equations 


£ 


3kn  (1  $kn)Xr 


jjb  1  nM 


Db  /v 


Ki  _  ph  (hi  -  Yk) 

— u k  v  ■ 

-&-n.  ^k 


where 


di  = 


EXn 

nb 


(14) 

(15) 


n^k  kn 


and  where  Dbkn  is  the  full-approximation  binary  diffusivity.  A  solution  for  can  be  obtained  through  an 
approximate  inversion16  as  follows 

where 


nM 

u kn 


■XkD%>, 


= 


1  +  Yk 
Xk 


DID 


D*k5kn  +  (1  -  Skn)-^  -  ( akD*k  +  anD*n )  + 

^ kn  1 


D*k  =  (l-Yk)Db, 
Mk. 

M 


&k 


=  ^(1+y‘)+Ey" 


D* 

^_n 

Db 

n^k  kn 


(16) 

(17) 

(18) 
(19) 


and  where  M  is  the  mixture’s  molar  mass.  To  avoid  in  the  above  method  artificial  singularities  for  mixtures 
with  vanishing  mass-fractions,  we  follow  the  procedure  of  the  EGLIB  multicomponent  transport  property 
library1'-18  .  To  this  end,  we  first  calculate  perturbed  mole  fractions  as 


XTr  =  Xn  +  e 


Yhn  Xn 
N 


~Xr 


(20) 


where  e  =  10  16  is  a  small  parameter.  Then  we  evaluate  the  perturbed  molar  mass  of  the  mixture  Mper  = 
XperMn  and  finally  the  perturbed  mass  fractions  as 


The  heat  flux  vector  reads 


\^per  _  Mn  -y'per 


q  —  — AVT  +  —  RuTapk)^~ 


(21) 


(22) 


where  A  is  the  thermal  conductivity,  Ru  the  universal  gas  constant  and  hk  the  partial  molar  enthalpy  of 
species  k. 
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Figure  1.  The  initial  potential  flow  entering  through  the  left  AB  and  right  A' B'  inlet  boundaries  and  exiting 
through  the  lower  AA!  and  BB'  boundaries.  The  solid  lines  represent  the  streamlines  and  the  arrows  the 
direction  of  the  flow. 


C.  Transport  properties 

Transport  properties  were  computed  according  to  the  most  up-to-date  methods14, 19-21 


D.  Equation  of  state 

All  NS  equations  are  coupled  with  an  EOS  which  is  selected  here  to  be  the  Peng-Robinson  (PR)  EOS22 


RUT 

p  = - r — 

VPR  Vmix 


& mix 

VPR  +  ^bmixVpR 


(23) 


from  which  T  and  p  are  obtained  as  an  iterative  solution  of  two  nonlinear  equations  that  satisfy  both  the 
values  of  p  and  e  as  obtained  from  the  solution  of  the  conservation  equations.  In  the  PR  EOS,  vpp  is  the 
PR  molar  volume  and  it  holds  v  =  Vpp  +  vs  where  vs  is  the  volume  shift  that  improves  the  accuracy  of  the 
PR  EOS  for  high-pressure  conditions23  ,  while  the  terms  amix  and  brnix  are  functions  of  T  and  Xj. 


III.  Initial  profile:  Inviscid  incompressible  potential  flow 

The  initial  flow  is  a  two-dimensional  potential  flow  in  which  two  opposing  streams,  each  of  constant 
density  and  composition,  are  injected  from  the  two  boundaries  AB  and  A'B'  of  length  Ly  in  the  ^-direction; 
the  streams  exit  the  domain  through  the  two  boundaries  AB  and  A'B'  of  length  Lx  at  the  y-direction  as 
illustrated  in  Fig.  (1).  At  the  point  C  of  the  left  inlet  boundary,  reference  inflow  variables  p'ref,T^epY-re^ 
and  plrej  are  specified  and  their  values  are  constant  with  time.  Similarly,  at  the  point  C'  of  the  right  inlet 
boundary,  reference  constant  in  time  inflow  variables  prrej ,  T'ej ,  Y[rej  and  prref  are  specified.  The  reference 
pressures  at  the  two  inlets  are  set  equal,  i.e.  plref  =  plref  =  pref,  so  that  for  the  initial  flow  described  in  the 
following,  the  stagnation  point  (S)  is  always  located  at  the  centerline  with  xs  =  Lx/ 2  and  also  ys  =  Ly/ 2. 

The  density  field  is  given  as 


P°(x,y) 


plref  if  0<x<Lx/2,\/y 
Pref  if  Lx/2  <  X  <  Lx,My 


(24) 


where  a  discontinuity  arises  to  p°(x,y)  at  x  =  Lx/2,Vy  when  plrej,  ^  Pref  -  The  species  initial  mass  fraction 
fields  are  given  as 


Y°{x,y) 


Yn,ref  lf  0<X<Lx/2,Wy 
Y^ref  if  Lx/2  <  x  <  Lx,My 


(25) 
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with  n  =  1  and  a  discontinuity  arises  to  Y®(x,y)  at  x  =  Lx/2,\/y  when  Y^  j,  ^  j..  The 

components  of  the  initial  velocity  field  are 


ul(x,y) 


2k(x  —  Lx/2) 


1  +  yPlef/Pref 
u°y(x,y)=  2<V~Lv/2)  , 

1  +  1 Jdef/Pref 

for  0  <  x  <  Lx/2,Vy ,  while  for  Lxj 2  <  x  <  Lx,My  are 

o  /  ni\ _  2k(x  Lx/2) 


ux(x,y) 


1  +  \J  Pre.fl  Pref 

2n(y  —  Ly/2) 


u"(x,y)=  v/ 


1  +  \J  Pref  /  Pief 


(26) 

(27) 


(28) 

(29) 


where  k  is  the  strain  rate.  ux(x,y)  is  continuous  V(x,  y),  while  uy(x,y)  has  a  discontinuity  at  x  =  Lx/2,\/y 
when  plref  y  prref. 

Using  (i)  the  common  reference  pressure,  pref,  of  points  C  and  C' ,  (ii)  the  Bernoulli  equation  which 
holds  for  inviscid  incompressible  flows  and  (iii)  the  fact  that  for  irrotational  (i.e.  potential)  flows  the  total 
pressure,  as  given  by  the  Bernoulli  equation,  is  the  same  for  all  points  of  the  flow  (otherwise  the  Bernoulli 
equation  would  only  hold  along  a  streamline),  we  obtain  the  initial  pressure  field,  p°,  as 


j  ^  r 

P°(x,y)  =  Pref  +  7fPlref  {u°x(xC,yc)  +  (u°y  (xC,  Vc))  ‘ 


rj  Pref 


{u°x{x,y))2  +  (uy(x,  y)Y 


(30) 


for  0  <  x  <  Lx/2,Wy,  where  xc  =  0  and  yc  =  Ly/2,  while  for  Lxj 2  <  x  <  Lx,\/y,  p°  is  given  by 


-L  "  2  r 

P°{x,y )  =Pref  +  ~ Pref  (u°x{xC' ,  VC')Y  +  (u°y  {xC> ,  VC’ ))  ' 


1  r 

2  Pref 


{u°x(x,y))2  +  (u0y{x,y)Y ]  ,  (31) 


where  xc  =  Lx  and  yc  =  Ly/2. 

By  construction,  the  initial  potential  flow  satisfies  the  steady-state  incompressible  Euler’s  equations  which 
in  “compressible"  form  can  be  expressed  as 


d(Ag)  ,  d(P°<)  n 

dx  dy 

d(p°u°xu°x)  d(p°u°xu°y)  dp0 

dx  dy  dx 

d(p°uxUy)  d{p°u°u°y)  dp0 

dx  dy  dy 

d{p°YSl)  ,  d(p°Yy°)  n  f_  1 

r\  I  0  —  vJ  [if'  —  J- 

ox  oy 


■■,N-  1) 


(32) 

(33) 

(34) 

(35) 


Vx,  y  except  for  the  locations  at  the  interface  line  of  the  two  streams,  i.e.  the  points  with  xs  =  Lx/2,\/y.  Equa¬ 
tions  (32)-(35)  will  be  used  next  as  a  basis  for  constructing  inflow  boundary  conditions  for  the  NS  equa¬ 
tions.  The  availability  of  p°(x,  y),  p°(x,  y),  Y®(x,  y), . . . ,  Y$(x,  y)  in  conjunction  with  the  EOS  allows  the 
determination  of  the  initial  temperature  field  of  the  flow,  T°(x,y),  by  solving 

p°{x,  y)  =  p(T°(x,  y),p°{x ,  y),Y°(x,  y), . . . ,  Y%(x,  y))  (36) 

for  T°(x,  y)  where  p°(x,  y)  is  given  by  Eqs.  (30)-(31),  p°(x,  y)  is  given  by  Eq.  (24)  and  Y®(x,y)  (n  =  1 . . . ,  N) 
by  Eq.  (25).  With  the  calculated  T°(x,  y),  the  initial  specific  internal  energy  field 

e°(x,y)  =  e(T°(x,y),p°(x,y),Y^(x,y), . . .  ,Y^(x,y))  (37) 

of  the  flow  is  calculated  as  well.  It  follows  that  the  initial  field  of  conservatives  variables,  U°(x,y),  can  then 
be  readily  constructed. 

Finally,  we  note  that  for  the  counterflow  configuration  the  use  of  a  potential  flow  as  initial  flow  is  consistent 
with  the  fact  that  even  when  starting  computations  with  a  plug  flow  (i.e.  zero  y-components  of  the  velocity 
vector  at  the  inlets),  the  flow  quickly  converges  to  a  potential  type. 
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IV.  Boundary  conditions 


A.  Inflow  boundary  conditions 

Instead  of  explicitly  imposing  timewise  constant  boundary  values  (also  known  as  hard  BCs)  for  the  conser¬ 
vative  variables  at  the  two  inlets  AB  and  A'B' ,  we  impose  the  steady-state  NS  compressible  conservation 
equations 
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ox  By 

■,N-  1) 

(41) 

for  p,  pux,  puy,  pY1, . . .  ,pYN_i .  Using  these  equations  as  BCs  means  that  the  initial  values  imposed  on  the 
inlets  AB  and  A'B'  for  these  variables  are  implicitly  constrained  to  be  constant  with  time.  Most  importantly, 
the  transverse  terms  that  give  rise  to  multi-directional  effects  at  the  inflow  boundaries  are  inherently  present 
in  the  BCs.  In  addition  to  the  above  BCs,  a  simple  algebraic  equation  of  the  form 


pet  =  pet(T(x,y),p{x,y),Y1(x,y), . . .  ,YN(x,y))  (42) 

is  used  for  the  boundary  values  of  pet  at  the  two  inlets,  where  et  is  calculated  using  the  current  values  of 
p,  pux,  puy,  pY\ , . . . ,  pYm-i  and  calculating  T  from  the  EOS  and  then  e,  so  that  the  values  of  pet  are  readily 
available. 

By  construction,  the  initial  values  of  p,  pux,  puy,  pY\, . . . ,  pYjv-i  at  the  inlets,  as  given  by  the  potential 
flow  of  Section  III,  satisfy  the  Euler  instead  of  the  NS  steady-state  equations.  Since  the  inlet  boundaries  are 
usually  located  away  from  the  flame,  these  initial  values  provide  a  good  initial  guess  for  the  NS  steady-state 
equations  which  can  be  corrected24  to  account  the  viscous  and  species  mass-flux  effects. 

B.  Outflow  boundary  conditions 

At  the  outflow  boundaries  A' A  and  BB' ,  the  solution  satisfies  the  unsteady  NS  equations  since  the  conditions 
there  are  the  outcome  of  the  processes  taking  place  in  the  interior  of  the  domain.  At  these  boundaries  we 
assume  that  p  stays  constant  with  time  as  initially  provided  by  the  Bernoulli  equation. 

V.  Numerical  scheme 

An  eight-order  explicit  finite-difference  scheme  is  used  for  the  spatial  derivatives.  After  spatially  dis¬ 
cretized,  the  steady-state  NS  equations  become  algebraic  equations  at  a  node,  whereas  the  unsteady  NS 
equations  become  ordinary  differential  equations  at  a  node.  As  a  result,  at  the  nodes  of  the  inflow  bound¬ 
aries  AB  and  A'B',  the  BCs  are  of  algebraic  type  since  Eqs.  (38)-(41)  and  Eq.  (42)  become  algebraic 
equations.  In  addition,  auxiliary  ghost  points  can  be  used  from  the  left  of  the  boundary  AB  and  the  right 
of  the  boundary  A'B'  (since  the  solution  is  known  there)  to  enhance  the  finite  difference  approximation  of 
the  spatial  derivatives  in  Eqs.  (38)-(41).  At  the  nodes  of  outflow  boundaries  A! A  and  BB' ,  the  BCs  are  of 
differential  type  and  one-sided  finite  differencing  is  used. 

The  resulting  numerical  system  of  differential-algebraic  equations25  is  integrated  in  time  with  the  differential- 
algebraic  solver  IDA  of  the  SUNDIALS  suite26 .  The  integration  method  used  in  IDA  is  a  variable-order, 
variable-coefficient  BDF  (Backward  Differentiation  Formula),  in  fixed-leading-coefflcient  form  where  the  or¬ 
der  of  the  method  varies  between  1  and  5.  The  BDF  method  can  handle  the  stiffness  introduced  in  the 
numerical  integration  of  the  NS  equations  due  to  the  presence  of  chemical  source  terms.  IDA  supports  paral¬ 
lel  computations  through  the  message  passing  interface  (MPI)  protocol  and  provides  a  routine  that  computes 
consistent  initial  conditions  from  a  users’  initial  guess24,2'  . 
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VI.  Reproducing  non-reacting  potential  flows  with  the  proposed 
differential-algebraic  boundary  conditions 

To  validate  the  proposed  differential-algebraic  BCs  of  Section  IV,  we  consider  two  symmetric  non-reacting 
potential  flows  where  either  H2  or  O2  enters  from  both  inlets.  The  goal  is  to  numerically  reproduce  these 
two  steady-state  potential  flows  at  computational  times  characteristic  to  reaction/diffusion  problems. 


Table  1.  The  inflow  conditions  at  the  reference  points  C  and  C'  of  the  left  and  right  inlet  (see  the  right  part 
of  Fig.  (1))  used  for  the  two  symmetric  Hi  potential  flows  cases  A  and  B. 


Case 

Species 

Pref  (kg/m3) 

Tref  (K) 

pref  (bar) 

«  (s-1) 

A 

h2 

0.08077 

300 

1 

2  x  103 

B 

h2 

0.08077 

300 

1 

4  x  103 

A.  The  H2  symmetric  potential  flow 

For  the  symmetric  H2  potential  flow,  two  situations  are  considered  with  different  strain  rates:  case  A  with 
k  =  2000  s-1  and  case  B  with  n  =  4000  s-1.  The  inflow  conditions  are  given  on  Table  1  and  the  size/meshing 
parameters  of  the  computational  domain  used  in  each  simulation  case  on  Table  2. 
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Figure  2.  Temporal  variation  of  p,vx,vy,p  and  Yh2  at  inlet  boundary  points  for  Case  A  of  the  symmetric  H2 
potential  flow  using  the  differential-algebraic  BCs  of  Section  IV.  The  inlet  conditions  are  shown  on  Table  1, 
the  size/meshing  parameters  of  the  computational  domain  on  Table  2  and  the  labeling  of  the  points  refers  to 
the  right  part  of  Fig.  (1). 


Regarding  case  A,  Fig.  2  shows  the  temporal  evolution  of  the  boundary  values  of  p,vx,vy  and  Yh2  at 
points  on  the  left  inlet  boundary.  The  prescribed  initial  boundary  values  are  excellently  maintained  by  the 
imposed  NS  steady-state  BCs,  Eqs.  (38)-(41),  at  the  two  inlets  AB  and  A'B' .  Moreover,  Fig.  4  shows  the 
initial  (t  =  0  ms)  spatial  variation  of  p.  vx,  vy.  p  and  Yh2  for  .T-sections  of  the  computational  domain  and 
the  spatial  variation  of  the  same  variables  obtained  at  t  =  16  ms.  Clearly,  the  steady-state  potential  flow  is 
accurately  reproduced  numerically  in  all  respects. 

Displayed  in  Figs.  3  and  5  are  the  equivalent  results  for  Case  B.  The  results  for  Case  B  are  equally 
excellent  in  that  the  potential  flow  is  maintained  even  for  the  larger  n  value. 
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Figure  3.  Temporal  variation  of  p,vx,vy,p  and  Yh2  inlet  boundary  points  for  Case  B  of  the  symmetric  H2 
potential  flow  using  the  differential-algebraic  BCs  of  Section  IV.  The  inlet  conditions  are  shown  on  Table  1, 
the  size/meshing  parameters  of  the  computational  domain  on  Table  2  and  the  labeling  of  the  points  refers  to 
the  right  part  of  Fig.  (1). 


Table  2.  The  lengths  Lx  and  Ly  of  the  computational  domain  (see  the  left  part  of  Fig.  (1))  and  the  number  of 
discretization  points  Nx  and  Ny  at  each  direction  used  for  the  simulations  of  the  two  symmetric  H2  potential 
flows  cases  A  and  B. 


Case 

Lx  [mm] 

Ly  [mm] 

N, 

Ny 

A 

10 

10 

104 

104 

B 

10 

10 

144 

144 

B.  The  O2  symmetric  potential  flow 

For  the  symmetric  O2  potential  flow  one  situation  is  considered  only.  The  inflow  conditions  are  given  on  Table 
3  and  the  size/meshing  parameters  of  the  computational  domain  used  in  the  simulation  on  Table  4.  The 
results  of  the  simulation  for  O2  are  illustrated  in  Figs.  6  and  7  which  represent  a  simulation  performed  with  a 
fluid  having  a  density  of  0(10)  larger  compared  to  that  of  H2  used  in  the  previous  two  simulations.  The  same 
high  fidelity  in  maintaining  numerically  the  imposed  initial  potential  flow  holds  as  in  the  H2  simulations. 

VII.  Summary  and  conclusions 

New  BCs  have  been  developed  for  the  counterflow  configuration  which  provide  a  very  concise  framework 
that  inherently  accounts  the  multi-directional  character  of  the  flow  at  the  boundaries.  Upon  discretization 


Table  3.  The  inflow  conditions  at  the  reference  points  C  and  C'  of  the  left  and  right  inlet  (see  the  right  part 
of  Fig.  (1))  used  for  the  symmetric  O2  potential  flow. 


Species  pref  (kg/m3)  Tref  (K)  pref  (bar)  k  (s  1) 
02  1.2837  300  1  2  x  103 
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Figure  4.  Spatial  variation  of  p,vx,vy,p  and  Yj. f2  for  Case  A  of  the  symmetric  H2  potential  flow  using  the 
differential-algebraic  BCs  of  Section  IV.  The  inlet  conditions  are  shown  on  Table  1,  the  size/meshing  parameters 
of  the  computational  domain  on  Table  2  and  the  labeling  of  the  sections  refers  to  the  right  part  of  Fig.  (1).  Left 
column:  f  =  0  ms.  Right  column:  t  =  16  ms. 
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Figure  5.  Spatial  variation  of  p,vx,vy,p  and  Yj. f2  for  Case  A  of  the  symmetric  H2  potential  flow  using  the 
differential-algebraic  BCs  of  Section  IV.  The  inlet  conditions  are  shown  on  Table  1,  the  size/meshing  parameters 
of  the  computational  domain  on  Table  2  and  the  labeling  of  the  sections  refers  to  the  right  part  of  Fig.  (1).  Left 
column:  f  =  0  ms.  Right  column:  t  =  6  ms. 
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Figure  6.  Temporal  variation  of  p,vx,vy,p  and  Yh2  inlet  boundary  points  for  the  symmetric  O2  potential  flow 
using  the  differential-algebraic  BCs  of  Section  IV.  The  inlet  conditions  are  shown  on  Table  3,  the  size/meshing 
parameters  of  the  computational  domain  on  Table  4  and  the  labeling  of  the  points  refers  to  the  right  part  of 
Fig.  (1). 


Table  4.  The  lengths  Lx  and  Ly  of  the  rectangular  computational  domain  (see  the  left  part  of  Fig.  (1))  and 
the  number  of  discretization  points  Nx  and  Ny  at  each  direction  used  for  the  simulation  of  the  symmetric  O2 
potential  flow. 


Lx  [mm] 

Ly  [mm] 

N* 

Ny 

10 

10 

200 

200 

of  the  steady-state  NS  equations,  the  BCs  at  the  inflow  boundaries  are  of  algebraic  type,  whereas  upon 
discretization  of  the  unsteady  NS  equations,  those  at  the  outflow  boundaries  are  of  differential  type.  This 
formulation  of  the  numerical  BCs  as  differential-algebraic  ones  requires  the  use  of  a  numerical  integration 
software  capable  of  handling  initial- value  problems  for  differential-algebraic  systems  of  equations  such  as  the 
IDA  of  the  SUNDIALS  suite26 . 

For  the  validation  of  the  differential-algebraic  BCs  two  symmetric  non-reacting  potential  flows  were 
considered  with  either  H2  or  O2  as  injected  species  from  both  inlets.  For  ED,  two  different  strain  rates  cases 
where  examined.  One  with  k  =  2000s-1  which  corresponded  to  an  injection  velocity  of  10  m/s  and  one 
with  k  =  4000  s-1  which  corresponds  to  an  injection  velocity  of  20  m/s.  For  O2,  a  single  case  was  examined 
with  k  =  2000  s-1  which  corresponds  to  an  injection  velocity  of  10  m/sec.  In  all  simulated  cases,  the  BCs 
were  able  to  maintain  accurately  the  imposed  boundary  values  at  the  inlet  boundaries.  Moreover,  the  initial 
potential  flow  was  reproduced  with  high  fidelity  over  the  entire  computational  domain  without  numerical 
artifacts.  For  the  aforementioned  cases  studied,  there  was  no  use  of  dissipative  filters  and  there  was  no  need 
of  introducing  relaxation  constants  into  the  BCs. 

Reactive  counterflow  simulations  with  the  differential-algebraic  BCs  are  currently  being  conducted  and 
will  be  the  subject  of  a  future  publication. 
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